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Liquid- vapor coexistence curves and critical parameters for hard-core 1:1 electrolyte models with 
diameter ratios A — a-/a+ = 1 to 5.7 have been studied by fine-discretization Monte Carlo methods. 
Normalizing via the length scale a± = |(cr+ + a-), relevant for the low densities in question, both 
T* (= kBTcO-±/q^) and p* (= PcO"±) decrease rapidly (from ~0.05 to 0.03 and 0.08 to 0.04, respectively) 
as A increases. These trends, which unequivocally contradict current theories, are closely mirrored by 
results for tightly tethered dipolar dimers (with T* lower by ~ 0-11% and p* greater by 37-12%). 



PACS numbers: 02.70.Lq, 05.70.Jk, 64.70. Fx 

The formation of coexisting fluid phases of different 
electrolyte concentrations in ionic solutions has been the 
topic of numerous recent experimental [Q , theoretical 
and simulation studies l^-^]- The simplest models for 
electrolytes — the "primitive models" — treat the solvent 
as a uniform dielectric continuum. The most studied sys- 
tem is the restricted primitive model (RPM) that consists 
of equisized hard spheres, half carrying a charge +q and 
half —q. Recent simulations |^-|^ agree with respect to the 
critical temperature and density for the vapor-liquid tran- 
sition, finding the remarkably low values T* ~ 0.05 and 
p* ~ 0.07 |§: see (|l|). By contrast to the RPM, the effects 
of charge and size asymmetry have not been extensively 
analyzed either theoretically or via simulation. 

Here we focus on the effects of size asymmetry on gas- 
liquid coexistence and critical parameters by studying 
hard-core primitive models for 1:1 electrolytes that have 
no restrictions on the relative magnitude of the diameters 
of the -|- and — ions, i.e., the so-called size-asymmetric 
primitive model (SAPM) The first claim to treat 

size asymmetry theoretically appears already in Debye 
and Hiickel's original paper |^,^ and extensions invoking 
Bjerrum ion-pairing have been analyzed Other mean 
potential approaches include the symmetrized Poisson- 
Boltzmann and modified Poisson-Boltzmann pQl s chemes. 
The mean spherical approximation (MSA) ]2|,|3|,p^ and 
hypernetted-chain (HNC) integral equations have also 
been applied. Currently, however, there are no simulation 
results available to check these various theories. 

This work, which extends provides a first study of 
the effects of size asymmetry on both critical parameters 
and liquid- vapor coexistence We find, in fact, a sys- 
tematic trend of Tc and pc with increasing size asymmetry 
that directly conflicts with the principal theories cited. 

To be specific, we consider a system of A'' hard spheres 
of diameter a+ carrying charges +q, and N of diame- 
ter (T_ carrying charges —q. The interaction energy be- 
tween two nonover lapping ions, i and j, of charges qi and 
qj (= ±g) separated by distance rjj is Uij — qiqj/Drij, 
where D represents the dielectric constant of the solvent 
which will be set to unity. The hard-sphere interactions 
are supposed additive so the (+,— )-ion collision diame- 



ter is (T± = ^ (cr+ + a-). This, in fact, provides the ba- 
sic length scale appropriate for defining both the reduced 
temperature and the reduced density via 



T* = kBTDa±/q^ and p* 



pal 



(1) 



1^,^ where p — 2N/V is the total ionic number density. 
Other definitions of the reduced density are, of course, vi- 
able and might be advantageous at high densities. 
However, at the low densities of interest here (2p* < 0.2) 
the formation of ion pairs, triples, chains and rings (see 
Fig. ^ below) is controlled almost exclusively by a±, which 
remains well defined even if or cr_ vanishes yielding 
point ions (idj. The reduced simulation box length is de- 
fined similarly via L* = L/a±. 

The key parameter for our study is the ratio of diameters 
of positive and negative ions, namely. 



A = cr_/(7-|_ 



(2) 



Because of symmetry with respect to the exchange of -|- 
and — ions, only A > 1 need be considered. 

In addition to the ionic systems, we have studied, for 
the sake of comparison , tightly tethered dipolar dimer 
systems consisting of N pairs of a positive and a nega- 
tive ion restricted to remain at separations ad satisfying 
Ci < CTd < 1.02a±. Interactions of the tethered dimers 
are otherwise identical to those of the ions. 

We adopt the methodology of |^2|. Neutral grand- 
canonical fine-discretization Monte Carlo simulations 
(characterized by a temperature, T, and a chemical po- 
tential for a pair of unlike ions, p) have been performed 
on cubic boxes of length L, under periodic boundary con- 
ditions. The positions available to each ion are the sites of 
a simple cubic lattice of spacing a. A "lattice refinement" 
parameter C, ^ a±la is introduced, so that when C, ^ oo 
the continuum is recovered. For values of C ^ 3, the RPM 
displays a gas-liquid coexistence curve that approaches the 
continuum case quite closely already for C = 5 ||l^. In a 
Lennard- Jones fluid studied using C = 10, the phase enve- 
lope and critical points of the lattice and continuum sys- 
tems were equal within the simulation uncertainties [ p^ . 
The structure of the liquid at short distances, as judged 
by the pair correlations, was also the same. 
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TABLE I. Critical parameters, T*[L*) and p*{L*), for the 
X = 1, — 10 ionic models with two values of eoo- (The la 
statistical uncertainties refer to the last decimal place.) 



TABLE IL Dependence on A of estimated critical parame- 
ters for (a) ionic and (b) tethered dimor models. 



L* 


lOOT*: = 1, 


eoo = oo. 


lOOp*: Eoo = 1, 


too = CO 


12 


5.09(1) 


4.97(1) 


7.0(3) 


8.0(2) 


15 


5.03(1) 


4.96(1) 


7.0(2) 


8.2(3) 


18 


5.00(1) 


4.96(1) 


7.4(2) 


7.9(2) 



We have thus adopted a refinement parameter oi( ~ 10. 
For the largest value of A we explore, namely A ~ 5|, the 
diameter of the smaller sphere is (T+ — 3a. When A = 1, 
typical correlation functions for both like and unlike ion 
pairs then agree to within graphical accuracy with the 
corresponding continuum results (for the same T*, p* and 
L*). For higher values of A (> 1), the discretization effects 
on the (+,— ) and (— ,— ) correlation functions decrease, 
while those on the (+, +) correlation functions increase, 
as CT-f/a becomes smaller; but the latter are very small at 
short distances because of the strong repulsions. 

The fine-lattice technique embodies precomputation 
and subsequent look-up of the Coulomb interaction be- 
tween any two lattice sites, including all periodic images. 
The Ewald sums were performed with conducting ("tin- 
foil") boundary conditions, i.e., Coo = oo; but for A =1, 
vacuum boundary conditions (eoo = 1) were also used to 
allow comparisons with ||]: see Table | and below. 

Biased insertions and deletions of pairs of unlike ions 
were performed for ionic models, following Q. Our teth- 
ered dimers, have 318 distinct configurations on the lat- 
tice. Dimers were inserted by randomly placing the — ion 
and selecting one of the 318 positions for the -I- ion. 

Histogram reweighting techniques were used to obtain 
the vapor-liquid envelopes up to T < 0.98Tc 0. Effective 
critical points for given L* were estimated using mixed- 
field finite-size scaling methods |l^ ], assuming Ising-type 
criticality. To discern a systematic dependence on A, this 
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FIG. 1. Phase diagrams for ionic systems with various 
size-asymmetries A — a-/a+. Statistical uncertainties are 
shown only when larger than symbol sizes. 
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approach should be satisfactory even though recent results 
[|9| (which indicate that the pressure should also enter the 
field-mixing) cast doubts on its full reliability. 

The two Ewald-sum boundary conditions for A = 1 
yield different critical values (see Table |) but extrapo- 
lation to L* = oo gives T* = 0.0495(2) for both cases and 
p* = 0.078(5) for e^o = 1 and 0.079(5) for e^o = oo. The 
agreement is excellent; but since the Coo = oo results vary 
less with L* they were used for the further simulations. 

These A = 1 results should approximate well those of 
the continuum RPM: recent studies yield 0.0488(2)- 
0.0490(3) for T* and 0.062(5)-0.080(5) for p*. Our 1% 
larger value of T* may be due to lattice discretization. 
Indeed, previous simulations [Q, found that both T*{L*) 
and p*{L*) estimates decreased slightly as C increased. 

The calculated phase diagrams for the ionic and teth- 
ered dimer systems are shown in Figs. |l| and ^ respec- 
tively. The critical points, which are listed in Table ||, 
were calculated Q using L* = 18 data while the subcrit- 
ical coexistence curves were obtained using L* — 12. The 
effects of size asymmetry are clearly strong, displaying a 
marked downward shift in T* and p* as A increases. (The 
lower values of T* result in smaller Monte Carlo accep- 
tance ratios and increasing sampling difficulties.) 

For A = 1, the critical temperatures of the ionic and 
dimer systems seem almost identical. Indeed, although 
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FIG. 2. Phase diagrams for tethered dipolar dimers. 
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FIG. 3. Reduced critical temperatures, , and den- 
sities, p*, as a function of the asymmetry parameter, 
a'(A) = (1 — A)^/(l + A^), as found by simulations (i) for ionic 
systems (open circles) and (ii) for tethered dimers (solid trian- 
gles); and as predicted by theory using (iii) the MSA (energy 
route) [2] (crosses) and (iv) extended Ebeling-Grigo theory (in 
the EG-Eb approximation) [3] (solid lines). 



Pc is about 37% higher, the overall phase behavior of the 
dimers is quite similar to that of the ionic systems, as 
stressed by Shelley and Patey When A increases, 

the critical temperatures of the dimers fall more rapidly 
but the critical densities approach those for ions, differ- 
ing by only 11% at A ~ 5.67: see Fig. ^ where these re- 
sults are depicted graphically vs the asymmetry param- 
eter w(A) = (1 — A)^/(l -I- A^). This parameter respects 
the symmetry under exchange of + and — ions and in- 
creases monotonically from — for the RPM, up to 
uj(oo) = 1, for point ions. Extrapolating our data to the 
point- ion limit, uj = 1, suggests a possibly common critical 
density, p*, for ions and dimers of ~0.015, with distinct 
critical temperatures, T*, in the ranges 0.020-0.025 and 
0.012-0.018. To check these speculations, however, may 
not be easy! 

Also shown in Fig. |[ as crosses, are the predictions of 
the MSA (using the energy route) for A = 1, 2, 10 and oo 
1^ . The absolute differences in T* and p* were to be antic- 
ipated (see, e.g., [^); but it is striking that the predicted 
changes with A in both T* and p* are opposite to those re- 
vealed by the simulations. The same failure to predict the 
correct trends is seen in the various calculations of Raincri 
et al. 1^ as illustrated by the solid curves in Fig. |^. These 
derive from extensions of the Ebeling and Grigo theory 
(which employs Bjerrum ion pairing). (See also ||20[] .) 



The conflict of our data with the theories cited is, per- 
haps, not so surprising when one recognizes the large de- 
gree of pair association that occurs already in the critical 
region of the RPM (A = 1) |2^-|2|. Indeed, the pres- 
ence of many such closely coupled dipolar pairs is what 
motivated the comparison with charged dumb-bells 
and our tethered dimer systems. But, as realized for some 
time |l^,^,^ and evident in sample configurations such 
as that in Fig. ^, dipolar systems undergo significant 
aggregation, primarily forming (+, — ) chains or "living 
polymers." Theories which mainly address the pair corre- 
lations cannot readily do justice to the geometrical aspects 
of the formation and interaction of such chains |2^ . 

Conversely, some insight into the lowering of T* as A 
increases may be gained by examining how the ground- 
state binding energies, = EbDa±/ ^ of various specific 
configurations depend on A. Thus for a neutral cluster of 
four ions (or two dipolar dimers) we find = 2.586 for a 
square "ring " when 1 < A < l-hA/2; but E^ falls smoothly 
when A exceeds 2.414 until, at A ~ 8.26, the lowest energy 
configuration switches from a planar diamond to a straight 



chain with E^ = 2.333. Similarly, if two long i 



chains 



are brought together, the excess binding energy per ion is 
El ~ 0.752 when 1 < A < 2.414 but decreases smoothly to 
0.698 when A grows larger. Other examples exhibit similar 
effects again rationalizing the observed drop in T*(A). 

The fact that tethered dimers show lower T* values (for 
A > 1) than ionic systems, supports the view that free ions 
assist in lowering the liquid free energy |Q. The result 
that the dimers require larger densities to stabilize the liq- 
uid is likewise consistent with this idea. In that connection 
the theory of pC| ] , in which dipolar dimers are solvated by 
-|- and — free ions (which strongly screen opposite ends 
of a dimer) may reasonably be regarded as approximating 
the solvation of a dipolar dimer by other dimers: these 
will screen by orienting in head-to-tail fashion. A better 
theory is much to be desired but, as various attempts il- 




FIG. 4. Snapshot of an ionic configuration with A = 5| and 
L* = 15 at T* = 0.0374 ~ imT^ and instantaneous density 
p* = 0.059 ~ 1.27/9*. 
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FIG. 5. Fraction, /, of ions in neutral clusters of n charged 
spheres in systems of size L* = 18 at p* ~ ip*. Solid symbols 
denote ionic systems, open symbols tethered dipolar dimers. 
The squares correspond to symmetric (A = 1) systems at 
T* = 0.050, the circles to A = 3 systems at T* = 0.045. 
From the top downwards (at n > 10) the reduced densities 
are p* = 0.026, 0.033, 0.019 and 0.027. 

lustrate (see, e.g., 0), that goal seems elusive. 

Finally, in an attempt to quantify some structural dif- 
ferences between ionic and tethered dimer systems, clus- 
ter densities were sampled for selected, comparable condi- 
tions: see Fig. |^. Gillan's definition of a cluster was 
used with a clustering distance Rq — Rc/cr± = 1.1. It is 
evident that there are more large clusters in the less sym- 
metric systems. But tethered dimers have much higher 
fractions of large clusters than do ionic systems, even al- 
lowing for the absence of charged clusters in the former. 

In summary, our fine-discretization simulations of hard- 
core 1:1 electrolyte models have provided unequivocal ev- 
idence that increasing the size asymmetry, measured by 
the diameter ratio A — a-jo-^.^ leads to sharp, mono- 
tonic drops in appropriately scaled critical temperatures 
and densities [see (|l|)]. Tightly tethered dipolar dimers 
(or dumb-bells ||l^) display broadly similar behavior but 
with relatively larger critical densities, and critical tem- 
peratures that decrease faster with increasing A. These 
trends are in severe disagreement with current theories 
and present what appear to be deep challenges to 
our theoretical understanding even at a qualitative level. 
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